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Abstract — We consider a certain class of large random ma- 
trices, composed of independent column vectors with zero mean 
and different covariance matrices, and derive asymptotically tight 
deterministic approximations of their moments. This random 
matrix model arises in several wireless communication systems 
of recent interest, such as distributed antenna systems or large 
antenna arrays. Computing the linear minimum mean square 
error (LMMSE) detector in such systems requires the inversion 
of a large covariance matrix which becomes prohibitively complex 
as the number of antennas and users grows. We apply the derived 
moment results to the design of a low-complexity polynomial 
expansion detector which approximates the matrix inverse by 
a matrix polynomial and study its asymptotic performance. 
Simulation results corroborate the analysis and evaluate the 
performance for finite system dimensions. 

I. Introduction 

Distributed antenna systems and large antenna arrays have 
recently attained significant research interest JT], 0. Both 
are considered as promising solutions to counter intercell 
interference and to increase the spectral efficiency of current 
cellular networks. Since these techniques rely in essence on 
a significant increase of the number of coordinated antennas, 
the computational complexity of the joint precoding/detection 
of the transmitted/received signals grows. This calls for low- 
complexity solutions. In this paper, we address this need by 
assessing the performance of a polynomial expansion detector 
adapted to the following general channel model. 

Consider a discrete-time N x K multiple-input multiple- 
output (MIMO) channel with output vector y 6 C N : 

y = Hx + n (1) 

where x = [x\, . . . , xk] T is the complex channel input vector 
satisfying E [xx H ] = l K , H = [hi • • • h K ] e C^*^ is the 
random channel matrix and n ~ CAf(0, o- 2 Ijy) is a vector of 
additive noise. The jth column hj £ C N of H is modeled as 

hj ■ = -^=RjW_j, j = l,...,K (2) 

where Rj € <C NxN is a deterministic matrix and the elements 
of Wj e C N are independent and identically distributed (i.i.d.) 
random variables with zero mean, unit variance and finite 
eighth moment. This channel model captures different types 
of wireless communication systems and generalizes several 
well-known channel models as discussed below: 



Distributed Antenna Systems: Let Rj = diag (fy, . . . ,TNj) 
with elements = y/Pj/d^J , where dij is the (normalized) 
distance between transmitter j and receive antenna i, /3 is the 
path loss exponent and pj is the transmit power of transmitter 
j. This model is suitable for distributed antenna systems HI 
where each transmitter sees a different path loss to each of the 
receive antennas since dij , . . . , d^j are different. 

Large-scale MIMO: Assume a receiver equipped with a very 
large antenna array (N 3> 1) as in [2|. Unless the antenna 
spacing is sufficiently large, it is likely that the received signals 
at different receive antennas are correlated. Our model allows 
to assign a different correlation matrix Rj to each transmitter. 

MIMO Multiple Access Channel (MAC): Consider a MIMO 
MAC from M transmitters equipped with K m , m = 1, . . . , M, 
antennas to a receiver with N antennas. Each point-to-point 
link has a different transmit and receive correlation matrix J4] : 

M 
m— 1 

where 3?r,i, . . . , <&r.m € C NxN are deterministic correlation 
matrices, * Ta e C* 1 ** 1 , . . . , * T ,M G £ KmxKm are 
nonnegative diagonal matrices, s C NxKl , . . . ,~W m G 
C Nx M are random channel matrices with i.i.d. entries with 
zero mean and variance 1/K, and Xi G C Kl , . . . , Xm G C Km 
are the transmit vectors. Let V _, K m = K. Setting Rj = 

ft* 3 e {l + ET=Ti 1 K l ,...,ET=iKi} and 
i = j — Y^iLi Kit we back to the model in d2l. 

In the sequel, we will study the asymptotic behavior of the 
moments /x„ of the matrix B = HH H , defined as 

/x„ = itrB", n = 0,1,2,... (3) 

under the assumption that N and K grow infinitely large at 
the same speed. In particular, we will derive deterministic 
approximations ~fi n of /x„ , such that /Lt„ — ~p n — > almost 
surely, for N, K — > oo. This result can be used, for example, to 
compute low-complexity approximations of the matrix inverse 
(B + cr 2 Ijv) _1 . The computation of this matrix arises in many 
practical applications, such as for linear multiuser detectors 
and beamforming strategies. We will focus exemplary on the 
linear minimum mean square error (LMMSE) detector. 



The LMMSE estimate x of x, assuming perfect knowledge 
of H at the receiver, is given as 



x = H H (B + a 2 I A r)- 1 y. 



(4) 



The computational complexity of this estimate is of order 
0(r 2 ) (6), where r = min(iV, K). A reduced complexity 
estimate can be obtained by approximating the matrix inverse 
in Q by the following matrix polynomial [ 3 ] 



(B 



<? 2 In) 



L-l 

E 

1=0 



(5) 



for some coefficients wi, where the filter rank L < r is chosen 
according to the allowable complexity. For a given transmitter 
k, the above polynomial expansion detector can be seen as 
a projection of y on the Lth Krylov subspace associated 
to the pair (B,hfc), i.e., the subspace of C N spanned by 
the vectors {h/^Bh/j, . . . ,B i_1 hfc}, and a weighting of the 
joint projections by the coefficients wi. Depending on L, 
the polynomial expansion detector achieves a performance 
between the matched filter (L — 1) and the LMMSE detector 
(L = r) [3 1 and allows, consequently, to trade-off performance 
for complexity. Moreover, <j3j allows for an efficient multistage 
implementation Q, Q, 0, where each stage I consists of a 
matched filter H H and subsequent "re-spreading" by the matrix 
H. In |8|, it was shown that the signal-to-interference-plus- 
noise ratio (SINR) at the filter output converges in certain 
cases exponentially in the filter rank L to the SINR output of 
the LMMSE detector. Thus, L does not need to scale with the 
system size to achieve close to optimal performance |9|. 

The optimal weight vector w = [wo • ■ • wl-i] T can be 
chosen to minimize the mean square error of the estimated 
vector x, i.e., 



w 



arg nun 

U=[«o,...,«L- 



E 



x-H 



L-l 



(6) 



The solution to this optimization problem is given as 

w = <fr~ 1 ip 
where $ € R^ xL and if 6 are defined as 

[<p]i = Mi- 

The computation of the weight vector w requires the calcu- 
lation of the moments fi%, . . . , which is still computational 
expensive for large L. However, under the assumption that the 
dimensions of H grow infinitely large, it was shown for several 
random matrix models (e.g. Q, (9j, iflOlO that the moments [i n 
can be closely approximated by their asymptotic counterparts 
p, n . These are independent of a particular realization of H 
and can be calculated based on the statistical properties of the 
channel matrix. If these properties change on a much slower 
timescale than the fast-fading channel fluctuations, the weight 
vector w can be precomputed using Jx n instead of /x n . Thus, 
the detector complexity depends only on the complexity of the 
projection on the Krylov subspace which is of order 0(r) J6|. 



Multistage or reduced-rank multiuser detectors were mainly 
considered in the context of code-division multiple-access 
(CDMA) systems as low-complexity solutions to the joint 
detection of a large number of user terminals with long 
spreading sequences Q. The asymptotic (universal) weight 
design was first studied in |7 1 for the equal transmit power case 
and then extended to more involved models, such as different 
transmit powers |9j, ifTTl . multi-path fading [ 1 1 and random 
unitary spreading sequences \ \2\. These results were then put 
on a common ground in |6| which compares different types 
of linear multistage detectors in terms of their complexity and 
asymptotic performance. Recently, also multistage detectors 
for asynchronous CDMA systems were considered in |[T3"1 . 

The asymptotic results in the above works are based on 
the almost sure (a.s.) convergence of the empirical spectral 
distribution (e.s.d.) of the matrix B to a compactly supported 
limit distribution. This limit distribution is in general given 
implicitly by its Stieltjes transform which can be computed 
based on the statistical properties of the underlying random 
matrix model. The asymptotic moments are then obtained by 
writing the Stieltjes transform as a moment generating function 
lfl4l Theorem 2.3] and relying on combinatorial arguments 
ifTOl or free probability theory |fl~2). 

The technique used in this work is different in two aspects. 
First, we do not require the existence of a limiting eigenvalue 
distribution of the matrix B. Instead, we provide for each pair 
(N, K) a deterministic approximation JL n of the moments \i n 
which becomes arbitrarily tight as N, K — i oo. Second, the 
moments are derived through iterated differentiation of the 
Stieltjes transform and can be computed by simple recursive 
equations. This is in contrast to ifTUl which requires an 
exhaustive search over complicated sets of indices. Hence, our 
results are more practical from an implementation perspective. 
Moreover, the asymptotic moments of the random matrix 
model |2} have not been considered in the literature before. 

The paper is structured as follows: Section [II] contains defi- 
nitions and related results. The asymptotic moments of B are 



(7) derived in Section III and the performance of the polynomial 



expansion receiver is studied in Section IV Numerical results 



are provided in Section [V] Section VI concludes the paper. 

II. Related results 

We need the following definitions and related results. De- 
note by "=^" and weak and almost sure convergence. 

Definition 1 (Empirical spectral distribution): Let A G 
^NxN ^ a Hermitian matrix with eigenvalues Ai, . . . , Xn. 
Denote F A the e.s.d. of A, defined as 



1 N 



Definition 2 (Stieltjes transform): Let F be a real measur- 
able function over M. with support Supp (F). For z € C \ 
Supp (F), the Stieltjes transform raf(z) of F is defined as 

f°° 1 

m F (z) = / dF(X). 

.1 _™ A — z 



Denote by S the class of functions / analytic over C\M + , such 
that, for z E C+, / G C+, zf E C+ and lirrij / _ i . 00 — iyf(iy) < 
oo. Such functions are known to be Stieltjes transforms of 
finite measures supported by IR + |14, Theorem 2.2]. 

Theorem 1 ( Kl~5\ Theorem 1]): Let D E <C NxN be a Her- 
mitian non-negative definite matrix and assume that D and 
the matrices Rj, j = 1,...,K, have uniformly bounded 
spectral norms (with respect to N). Let N, K — > oo, such that 
< liminf |C < limsup |? < oo. Then, for any z E C \ R+, 



^trDCB-^) 



- ^trDT(z) 



where T(z) g C JVxiV is defined as 



T(z) 



K 



if ^ 1 



- zl 



(9) 



and the following set of if implicit equations 
1 



K 



trRjR^T(z), j = l,. 



admits a unique solution (Si(z), . . . , 5k (z)) E S K . Moreover, 
denote by F the distribution function whose Stieltjes transform 
is given by m(z) = -^trT(z). Then, almost surely, 

,F B - F => 0. 

III. Asymptotic Moments 

In this section, we state our main results. The proofs of 
Theorems [2] and [3] are provided in the appendix. 

Theorem 2: Let F be the distribution function as defined in 
Theorem [T] and denote by JI , JL 11 . . . the successive moments 
of F, defined as ~p n = X n dF(X). These moments can be 

(-1)" 1 



calculated as 



M„ = 



;trT„ 



n\ N 

where T„ is defined recursively by the following set of 
equations for n > 0: 



-n+l — 



EE 

i=0 j=0 



i Qi 



Qr 



K 



1 K 



,RfcR 



H 



fc=l 



fk, 



n+l 



where T 



=EE 

= -^trR fe R£T 

: IjV, fkfl 



(n-i + l)fk,jfk,i-jh,' 



-1 and 4 n 



itrR fc R£ Vfc. 



Remark 3.1: While Theorem [2] allows to compute the mo- 
ments Jl n of F, it does not imply the a.s. convergence of 
fi n and ~p n in general. Theorem [3] provides some sufficient 
conditions for which this convergence holds. 



Remark 3.2: Although difficult to show analytically, one 
can verify numerically that Theorem [2] coincides with [10 
Theorem 1] for Rj = diag(rij, . . . , r/vj), j = 1, . . . , K , 

If the matrices Rj are drawn from a finite set of matrices, 
we get the following stronger result: 

Theorem 3: For fixed M > 0, let 1Z = {Rl, . . . , Rm} be 
a set of complex N x N matrices and let D E <C NxN be a 
non-negative definite Hermitian matrix. Assume that D and 
R m , m = 1 , . . . , M, have uniformly bounded spectral norms 
(with respect to N). Let R 3 E 1Z for j — 1,...,K, Assume 
N,K — > oo, such that < liminf^ < limsup^ < oo. 
Then, for n = 0, 1, 2, ... , 



1 (-1)" 1 
trDB™ ^ trDTV 







N n! N 

where T n is given by Theorem [2] This implies in particular, 



Mr, 



0. 



Loosely speaking, Theorem [Tj states that, for large matrix 
dimensions, the e.s.d. F s of the matrix B can be closely 
approximated by a deterministic distribution function F. Thus, 
the optimal weighting vector w can be approximated by 
replacing the moments \i n of F B in ([H} by the moments ~p n 
of F. Using the result of Theorem [2] we can compute an 
approximate weight vector w = [wq . . . wl-i] as 



w = <& 1 <p 

where $ E M.+ x L and Tp E are defined by 



(10) 



[*]« =Vi+j +° 2 [h+j-i (11) 
W\i =¥%■ 

IV. Asymptotic Performance Analysis 

We consider now the asymptotic performance of the polyno- 
mial expansion receiver in terms of the received SINR 7/. for 
a given transmitter k. With weight vector w, the fcth element 
£k of the estimated vector x reads 



L-l 



Xk 



h£ w i Bl ( Hx + n ) 



(12) 



1=0 



One can easily show that the associated SINR 7^ can be 
expressed as |6l Eq. (18)] 



Ik 



where $fc G M^ xi and (p k E are given as 



[*fc], 



[B I+J ] 



fcfc 



[B 



i+j-ii 



1 fcfc 



(13) 



(14) 



[fk]i = [B l ] 



fcfc 



The next theorem provides a tight deterministic approximation 
of the terms [B"] fefc = h^B n_1 hfc in the asymptotic limit. 



Theorem 4: Under the assumptions of Theorem [3] the fol- 
lowing convergence holds: 



\B n ' 



k-k- 







where 



i=0 



Mn - y^Mn-i-i " Tj~ — trR fc RfcTj 



n > 1 



and T n is given by Theorem [2] The initial values of the 
recursion are JIq = 1 and To = I at. 

Proof of Theorem K The proof follows the same steps 
as @ Theorem 1] and will not be given here. ■ 

Replacing [B"] fe/c in ( fl4] > by ~pf^ and w in ([13) by w, we 
can obtain a deterministic approximation of the SINR 7^ at 
the output of the polynomial expansion receiver. 

V. Numerical Results 

Consider a MAC from K = 40 single-antenna transmitters 
to a receiver with N = 100 antennas. We use an extended 
version of Jake's model [4] for the generation of the matrices 



Rj-. Let Rj = & J 



1/2 



and &j € 



iNxN 



-"max 



9 ■ 

mm 



exp 



be defined as 
27ri 

T 



-dki cos(x) ) dx 



where dki — 2A(fc — Z) and (f> 3 mia , ^max are drawn independently 
from the intervals [— ir, 0] and [0, ir], respectively. The interval 
I^min' ^max] can be seen as the angular spread of the signal 
from transmitter j, X is the wave length, and dki is the spacing 
between the receive antennas k and I. We assume Rayleigh 
fading channels, i.e., in |2]) are independent standard 
complex Gaussian vectors. The covariance matrices &j are 
chosen at random at the beginning and then kept fixed while 
we average over many realizations of the channel matrix H. 
We denote by SNR = l/c 2 the transmit signal-to-noise ratio. 

Fig.[T|shows the average received SINR £[7^.] of a randomly 
chosen transmitter as a function of the SNR for the matched 
filter, the LMMSE detector and the polynomial expansion 
detector with approximate weights for L = {2, 3, 6}. Markers 
correspond to simulation results and solid lines to the deter- 
ministic SINR approximations. The error bars indicate one 
standard deviation of 7/c in each direction. Similar to [16], the 
asymptotic SINR of transmitter k for the LMMSE detector 
can be easily shown to satisfy 

1 



—lmmse = trR fc RfeT(-l/SNR) 
K 



where T(z) is given by Theorem [T] We observe a good fit 
between the deterministic approximations and the simulation 
results for the average SINR. However, the standard deviation 
of the SINR increases with L. This is because the higher order 
moments converge slower to their deterministic approxima- 
tions and exhibit therefore stronger fluctuations. Nevertheless, 
the average SINR performance of the polynomial expansion 
detector with L = 6 is already close to the performance of the 
LMMSE detector. 
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Fig. 1 . Average received SINR versus SNR at the output of the matched filter, 
LMMSE detector and the polynomial expansion detector with approximate 
weights for different values of L. Markers correspond to simulation results, 
solid lines to the deterministic SINR approximations. Error bars indicate one 
standard deviation in each direction. 
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Fig. 2. Average theoretical bit error rate versus SNR for the matched filter, 
LMMSE detector and the polynomial expansion detector with approximate 
weights for different values of L. 



Fig. [2] depicts the theoretical average bit error rate (BER) 
over SNR for the different detectors. Assuming binary phase- 
shift keying (BPSK) modulation and Gaussian interference, 
the BER is given as E [(5(^/7^)] where Q(x) is the Gaussian 
tail function. We can clearly see a performance increase of 
the polynomial expansion detector with L, although the BER 
saturates at high SNR. Although not explicitly shown here, 
one can even observe a performance decrease for large values 
of L. As mentioned before, this is due to the low accuracy of 
the approximate weights caused by a slow convergence of the 
higher-order moments to their deterministic approximations. 



VI. Conclusion 

We have derived asymptotically tight deterministic approx- 
imations of the moments of a certain class of large random 
matrices, useful for the study of distributed antenna systems 
and large antenna arrays. We have applied these moment 
results to the design of a polynomial expansion detector which 
significantly reduces the computational complexity of mul- 
tiuser detection compared to the LMMSE detector. Moreover, 
we have derived an explicit expression of the asymptotic SINR 
at the output of this detector and verified its accuracy and 
performance for finite system dimensions by simulations. 

Appendix 

Proof of Theorem |2f From Definition [2j it is easy to 
see that the moments fi n of the distribution function F can 
be obtained through successive differentiation of the function 
~m(—~), i.e., 



(-1)™ d n 
n! dz n 
(-1)" d n 
n! dz n 

(-1)™ r d n 

dz r 



1 ( 1 

— m — 

z \ z 



1 

z\+ 1 
1 



z\ + 1 



z=0 

dF(X) 
dF(X) 



z=0 



z=0 



X n dF(X) 



where we could exchange the order of differentiation and 
integration since l/(zX + l)dF(X) is analytic for z > 0. 
Consider now the following function for z > 0: 

rj(z) = -to [ — 
z \ z 

and denote rj n (z) its nth derivative with respect to z. From 
Theorem [T] we have 

f \ 1 ( 1 
riiz) = —m — 

z \ z 



N { k U 1 + s A-j) 

1/1 RjR 1 ^ 



+ z5j,o(«) 



1 

iV 



trTo(z) 



where 



and ((5 10 (z), . . . , Jk-.o(z)) € is the unique solution to the 
K implicit equations: 

S jfi {z) = itrR^ToOO, j = l,...,K. 



Denoting T„(z) = J » » we have 
= ^trT„(z). 

In order to find the derivatives T„(z), we need the following 
additional definitions. For k € {1, ... , if}, let 



h,o{z) 



1 



1 + .9fc,o( z ) 
ifc,o( z ) = zfkfl(z) 

and denote 4,n(z), 9k,n{z), fk, n (z), and t fe ,„(» their nth 
derivatives, respectively. Furthermore, let 

1 ^ 

= ^H^,oWRfcRfc 

fe=i 

and denote Q n (z) = d • We continue by writing 

Ti(z) =To(g) Qi(^)To(z) . (15) 

= G (z) 

From the Leibniz-rule for the nth derivative of the product of 
two functions^ we have 



T n+1 (z) =J2 ( ■ )T n -i(z)Gi(z), n > 

GhW^Wq^iWIz), n>0 
where G„(z) = d ■ Replacing the last equation in the 



second last yields 

n i 



i=0 j=0 



T„- l (z)Q l _, +1 (z)T,(z). 



Straight-forward differentiation of Qo(z) leads to 
1 - 

Qn(«) = — 5^t fc)Tl (a:)RfcRjJ, n > 0. 



(16) 



(17) 



k=l 



The last step is to find explicit expressions of tk, n (z)- From 
the Leibniz-rule, we have 

tk,n(z) = nf k , n -i{z) + zf k ,n{z) , n > 0. 
Consider now fk.i( z ) the first derivative of fk,o{z): 

(1 + g k>0 ) v — - 
The higher order derivatives are calculated as 



fk,n+l(z) = ^ \ n )rkA z )9k,n-i+l(z) 



'For two functions u(x) and v(x), 

n ln\ d n ~ l u(x) d 1 v(x) 



d n (u(x)v(x)) 
dx n 



En rn\ d u{ 
i=0 \i) 



dx 1 - 



where 



rk. 



i=0 



Combining the last two equations yields 



fk,n+l(z) = 



fk.j (z)fk,i-j (z)gk,n-i+l(z) 



(18) 



where gk, n (z) can be easily calculated as 

9k,n( z ) = n $k,n-l{z) + Z6 k>n (z) 

and Sk, n {z) is given by 

tik,- 



= 4trR fe RfeT„(z). 

A 



(19) 



Since we are only interested in the case z = 0, we will drop 
from now on the dependence on z and write, e.g., T n instead 
of T„(0). In this case, the expressions of Qk,n{z) and tk, n (z) 
simplify to 



9k,n 
^k,n 



nSk, n -l 
n fk,n-l- 



Replacing these quantities in ( fT7| ) and ( [IS) , together with ( fT6| ) 
and ( fT9) i leads to the desired result. Note that T — In, 
fk.o = — 1 an d $k,o = ^trRfeR-fc- Moreover, T n+ i depends 
on T , . . . , T„ and Qi, . . . , Q n +i. Since Q n +i depends only 
on f k<0 , f k<n and / fe) „, T n+i can be recursively calculated 
from the given initial values. ■ 

Proof of Theorem ^ Both J?trD (B — zlpj)^ 1 and 
-^trDT(z) as defined in Theorem fll are Stieltjes transforms 
of finite measures which we denote by ir and W, respectively. 
Thus, Theorem [T] also implies that, almost surely, 

7T — 7f => 0. 

Similar to the proof of Theorem[2]we can express the moments 
of 7r and 7f as 



\ n ir{d\) 



(-1)" d" 



1 1 



i! dz n \ z N 



trD B 
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-trDB 1 



and 



A n 7f(dA) 



(-1)" d" /l 1 
(-1)" 1 



trDT(-l/z) 



2 = 



i! TV 



tiDT n . 



The support of n is almost surely compact as D has bounded 
spectral norm and the spectral norm of B is almost surely 



bounded due to the following inequalities: 

M 



|B|| < J2 W*mK 
1 



rn—1 

< MR sup 



W W 



A' 

A". 



W W H 

vv m VT m 



A/ A sup I 1 



for some R > sup„ 



N 



and W m g ^ e j n g ran( j om matrices with i.i.d. ele- 

ments with zero mean, unit variance and finite eighth moment. 
The almost sure convergence of the spectral norm in the last 
step follows from ifTTIl . The almost sure weak convergence of 
7r and 7f implies by Theorem 25.8 (ii)], that 



< oo 



/(A)Tr(dA) - / /(A)Tf(dA) ^ 



(20) 



for any bounded, continuous function. Since the support of ir is 
almost surely bounded and the support of n can be shown to be 
bounded following similar steps as in |4| Proof of Theorem 2, 
Part B], the convergence in ( |20| > also holds for any continuous 
function. Choosing /(A) = A™ concludes the proof. ■ 
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